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ABSTRACT 

We have conducted three-dimensional self-gravitating radiation hydrodynamical models of 
gas accretion onto high mass cores (15-33 Me) over hundreds of orbits. Of these models, one 
case accretes more than a third of a Jupiter mass of gas, before eventually undergoing a hy- 
drodynamic collapse. This collapse causes the density near the core to increase by more than 
an order of magnitude, and the outer envelope to evolve into a circumplanetary disc. A small 
reduction in the mass within the Hill radius (Ru) accompanies this collapse as a shock prop- 
agates outwards. This collapse leads to a new hydrostatic equilibrium for the protoplanetary 
envelope, at which point 97 per cent of the mass contained within the Hill radius is within the 
inner 0.03 Rh which had previously contained less than 40 per cent. Following this collapse 
the protoplanet resumes accretion at its prior rate. The net flow of mass towards this dense 
protoplanet is predominantly from high latitudes, whilst at the outer edge of the circumplan- 
etary disc there is net outflow of gas along the midplane. We also find a turnover of gas deep 
within the bound envelope that may be caused by the establishment of convection cells. 

Key words: planets and satellites: formation - methods: numerical - hydrodynamics - radia- 
tive transfer 



1 INTRODUCTION 

The ideas discussed in this paper begin with |Perri & Cameron] 
{1974}, who stated that "when the mass of the core becomes suf- 
ficiently great, the surrounding gaseous envelope will become hy- 
drodynamically unstable against collapse onto the planetary core". 
This process is controlled by the battle between gravity that acts to 
contract an envelope onto the core and the gas pressure which acts 
to support it. 

Mizuno ( 19801 performed stability calculations to determine 
the combinations of core masses and opacities that would make 
a protoplanetary envelope unstable to collapse. Work of his con- 
temporaries considering the structure of giant planets in the solar 
system suggested that each such planet (Jupiter, Saturn, Uranus, 
& Neptune) possessed a solid core with a mass of order 10 M e 
l |Slattery|1977"l|Hubbard & Macfarlane|1980) . Using these values 
as a target, and assuming a fixed accretion rate, |Mizuno| concluded 
that a grain opacity of K « 1 cm 2 g _1 was required in the enve- 
lope material during formation to trigger a collapse when the core 
mass was « 10 M ffi . Lower opacities were found to lead to enve- 
lope collapse at lower core masses. Following the envelope col- 
lapse [Mizuno] states that continued accretion is likely required for 
the protoplanet's to attain their final masses. 

It was suggested later by Bodenheimer & Pollack ( 1986 1, who 
performed evolution calculations that included evolution beyond 
the attainment of the critical core mass, that rather than a dynam- 
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ical collapse, the envelope may quasi-statically contract onto the 
solid core. They suggested that if a sufficient mass of molecular 
hydrogen in the envelope was apt to undergo dissociation, then this 
would remove enough energy from the contraction to bring about a 
dynamical collapse. However, their models indicate that such disso- 
ciative regions possess an insufficient fraction of the envelope mass 
for this to occur. 

In the early 1990s Wuchterl wrote a series of papers exploring 
the evolution of a protoplanetary envelope through the hydrostatic 
phases approaching the critical core mass, and in what he found to 
be a subsequent hydrodynamic phase (Wuchterl 1991 1. Solving the 
equations of radiation hydrodynamics in one-dimension, Wuchterl 
indicated that following a period of quasi-static contraction, during 
which the envelope heats up, the transport of this heat out through 
the envelope by convection and radiation perturbs the hydrostatic 
equilibrium. In particular Wuchterl cites the k - mechanism as the 
means of exciting the dynamical waves that destabilise the enve- 
lope. The result of this hydrodynamic evolution was the ejection of 
a large fraction of the envelope mass, rather than an inwards col- 
lapse as had previously been supposed by |Perri & Cameron ( 1974| >. 

|Pollack et al.| ( |1996| > continued performing models assuming a 
quasi-hydrostatic contraction, and suggested that further work was 
required to consider the hydrodynamic evolution of young proto- 
planets to establish the proper evolution scenario. A step in this di- 
rection was taken by Tajima & Nakagawa (1997) who performed a 
stability analysis of a growing protoplanetary envelope using a dis- 
tinct numerical code to that of previous authors. Their aim was to 
determine whether quasi-static contraction, or an envelope instabil- 
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ity akin to that suggested by Wuchterl was the more likely evolu- 
tionary course for a growing protoplanet. They perturbed the enve- 
lope at intervals during its evolution to see if such action might push 
a marginally stable system towards instability. They concluded that 
quasi-static contraction was viable for a protoplanet growing all the 
way to a Jupiter mass. 

There has since been a substantial amount of work consid- 
ering the gas accretion rates that cores might achieve in circum- 
stellar discs with a variety of properties. |Ikoma et al.| {2000} per- 
formed quasi-static evolutionary models to determine the depen- 
dence of gas accretion upon core mass, grain opacity, and the core's 
accretion history, finding that these factors were strongly, mod- 
erately, and weakly significant respectively. |Bryden et al.| J 1999} 
and [Lubow et alT] ( |1999| > performed locally-isothermal hydrody- 
namics simulations of discs containing planetary cores, the latter 
finding that the accretion rate drops off as the protoplanet mass be- 
comes very large; a result of the broadening disc gap that it forms. 
Further hydrodynamic models with more realistic thermodynamics 
were performed by|rPAngelo et al.|12003a},|Klahr & Kley|((2006) , 
|Paardekooper & Mellema| ( |2008| l, and |Ayliffe & Bate| ( |2009a) , fiiid- 
ing similar turn overs in the accretion rate with increasing mass, and 
illustrating the impact of grain opacity on accretion. These models 
have generally had limited resolution in the vicinity of the proto- 
planet, and though Ayliffe & Bate (2009a I achieved high resolu- 
tion, the evolutionary period was extremely short due to the com- 
putational demands of the calculations. 

In this paper we report results from three-dimensional self- 
gravitating radiation hydrodynamics calculations that resolve the 
protoplanet's envelope, whilst modelling its hydrodynamic evolu- 
tion and growth within a section of a circumstellar disc. Using high 
mass discs (though still stable; Toomre Q » 1), and assuming low 
opacities, we achieve accretion rates that allow significant envelope 
growth in only a few hundred orbits, allowing us to examine their 
development. Our computational method is described in Section|2] 
followed by our results in Section[3] and a discussion of their rela- 
tionship with previous results in Section]?] 



2.1 Model setup 

The calculations we have performed were conducted by modelling 
a small section of a circumstellar disc, centred upon a protoplane- 
tary core and corotating with its orbit. The protoplanet is modelled 
as a gravitating mass with a 'surface'. It attracts gas from the disc, 
building up an atmosphere which is supported by the surface. The 
disc section measures r = 1±0. 15 r p (5.2±0.78 AU), and0 = ± 0. 15 
radians. The protoplanet is sited at a radius of r p , which in all cases 
is equivalent to 5.2 AU, and orbits about a star of 1 M Q . The disc 
has a surface density profile of £ oc r~ l/2 , and a temperature profile 
of T g oc r (giving a scaleheight of H/r = 0.05), where these pro- 
files are equivalent to those used in our previous work, and were 
originally chosen to match work of |Lubow~ et al. (1999) and |Bate| 
|et al.| ((2003 1. The initial temperature at r p is == 73 K, where this tem- 
perature is taken to principally result from stellar irradiation. Near 
the midplane of the disc, this initial temperature may increase due 
to viscous heating as the model evolves if the opacity is sufficient 
to insulate the region against rapid radiative cooling. However, at 
the upper and lower radiation boundaries (see section [23] below), 
where the medium is always optically thin, an assumption of tem- 
peratures dominated by the stellar irradiaton is reasonable. The sur- 
face density at r p has an unperturbed value of 750 g cirr 2 , which 
gives a relatively massive disc, comprising 0.1 M Q of gas within 
25 AU. We use a high disc surface density because Ayliffe & Bate 
(2009a) found that increasing the disc surface density resulted in 
somewhat faster accretion rates onto embedded protoplanets and 
the goal of this paper is to investigate the three-dimensional evolu- 
tion of protoplanets that accrete massive gaseous envelopes. 

Further details of our model setup are given below, but the 
method is identical to that employed in Ayliffe & Bate (2009a), 
and that paper contains somewhat more extensive information, in- 
cluding resolution tests. We note that the most interesting model 
discussed in this paper, a case which results in a hydrodynamic 
envelope collapse, required more than 6 CPU years to reach its fi- 
nal state. It is this extensive calculation time which has limited the 
number of models that is has been possible to perform. 



2.2 Disc sections 



2 COMPUTATIONAL METHOD 

The calculations discussed in this paper were performed using a 
three-dimensional SPH code. This SPH code is derived from a code 
first developed by |Benz| ([T990 ; Benz et al. 1990) which has un- 
dergone substantial modification in subsequent years. Energy and 
entropy are conserved to timestepping accuracy by use of the vari- 
able smoothing length formalism of Springel & Hernquist ( 2002 1 
and Monaghan (2002), where our particular implementation is de- 
scribed in |Price & Monag han"] j2007| >. Gravitational forces are cal- 
culated and particle neighbours are found using a binary tree. Ra- 
diative transfer is modelled using the flux-limited diffusion approx- 
imation, employing the method developed by |Whitehouse, Bate, &| 
|Monaghan| \2005) and |Whitehouse & Bate] l |2006} . Integration of 
the SPH equations is achieved using a second-order Runge-Kutta- 
Fehlberg integrator with particles having individual timesteps {Bate] 
|1995[ >. Gas within the models is subject to an artificial viscosity, im- 
plemented in a parameterised form as developed for SPH by |Mon-| 
|aghan & Gingold] ( 1983}, and modified to deal with high Mach 
number shocks by Monaghan] ljT992;). The code has been paral- 
lelised by M. Bate using OpenMP and MPI. 



Within the disc section being modelled, the particles are initially 
distributed according to the underlying density profile, and their 
velocities are Keplerian. The protoplanet is orbiting its star in an 
anticlockwise direction, and in the corotating frame of the mod- 
elled section this leads particles at r < r p to orbit anticlockwise, and 
those at r > r p to orbit in a clockwise fashion. Particles encounter- 
ing the boundary of the section are removed from the calculation, 
whilst a number of ghost particles beyond the domain of the calcu- 
lation, along its boundaries, act to replicate the pressure and viscous 
forces expected from a continuous disc. To prevent the depletion of 
gas within the disc section, particles are injected along the bound- 
aries where the material should flow in. For the anticlockwise orbit- 
ing gas this input is along the <§> = -0.15 radian boundary, between 
r = 0.85 - 1.0 r p , whilst the clockwise flowing gas is injected along 
the (f> = 0.15 radian boundary, between r = 1.0 - 1.15 r p . The in- 
jection scheme does more than simply replace gas which leaves the 
section. The velocity and density structure of the injected gas is ob- 
tained from three-dimensional global simulations of protoplanets 
embedded in discs performed using ZEUS (Bat e et al.|2003) , such 
that a gap is opened in the disc which corresponds to the mass of 
the embedded protoplanet. The models of |Bate et al.j used a sur- 
face density of 75 g cirr 2 , but the particles injected in the calcula- 
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tions discussed here have their masses scaled to reflect the chosen 
surface density of 750 g cnT 2 . As the protoplanet grows the gap 
width is suitably increased by interpolation through the protoplanet 
mass range provided by |Bate et al.H2003) (1 M e -1 Mj upiter ). The 
injected particles come to dominate the section after less than five 
orbits, ensuring the structure is consistent with the presence of the 
protoplanet. 



Model Core mass Core radius X Physical Opacity 





[M ] 


[Mffi] 


[i'p] 


core size 


[% IGO] 


A 


4.54 x 1(T 5 


15 


2.2 x 10~ 4 


10 


0.1 


B 


6 x ur 5 


20 


6.6 x 10~ 5 


3 


0.1 


C 


6 x ur 5 


20 


6.6 x 10~ 5 


3 


1 


D 


6 x ur 5 


20 


2.2 x 10~ 4 


10 


0.1 


E 


6 x i(r 5 


20 


2.2 x 10~ 4 


10 


1 


J 


l x ur 4 


33 


2.54 x 10~ 4 


10 


0.1 



2.3 Thermodynamics 

All the calculations discussed in this paper were performed using 
radiation hydrodynamics, employing a two temperature (gas and 
radiation) radiative transfer scheme using a flux-limited diffusion 
approximation (as described by|Whitehouse et al. 2005 and White- 
|house & Bate 2006). Work and artificial viscosity act to increase 
the thermal energy of the gas, and work done on the radiation field 
increases the radiative energy, which can be transported via flux- 
limited diffusion. The energy transfer between the gas and radiation 
fields is dependent upon their relative temperatures, the gas density, 
and the gas opacity. 

The gas is treated using an ideal gas equation of state p = 
pTgR g /fi where R g is the gas constant, p is the density, T g is the gas 
temperature, and p is the mean molecular mass. The thermal evolu- 
tion takes into account the translational, rotational, and vibrational 
degrees of freedom of molecular hydrogen (assuming a 3:1 mix of 
ortho- and para-hydrogen; see Boley e t al.|2007) . Also included are 
molecular hydrogen dissociation, and the ionisations of hydrogen 
and helium. The hydrogen and helium mass fractions are X = 0.70 
and ¥ = 0.28, respectively, whilst the contribution of metals to the 
equation of state and the thermal evolution is neglected. 

The flux-limited diffusion scheme transfers energy between 
SPH particles, which does not enable it to radiate into a vacuum. 
In order that the disc can cool from its upper and lower surfaces, 
a boundary is applied that maintains the initial temperature profile 
in the high atmosphere of the disc. This boundary is situated at a 
height above/below the midplane that corresponds to the edge of 
the optically thick region, that is where the optical depth (t) from 
outside the disc to that depth is r as 1. SPH particles comprising 
the boundary regions evolve normally, but their energies are set 
according to the initial radial profile, allowing them to act as energy 
sinks. 



2.4 Opacity treatment 



We use the opacity tables of Pollack et al. < |1985|> to provide grain 
opacities, whilst the tables of Alex ander] J 1975) (the IVa King 
model) provide the gas opacities at higher temperatures when the 
grains have sublimated. The former table gives interstellar grain 
opacities (IGO) for solar metallicity molecular gas, but in this work 
we reduce these opacities by orders of magnitude below this nomi- 
nal level; we divide by factors of 100 and 1000. The justification for 
this comes of the likely agglomeration or sublimation of grains in 
the vicinity of a forming protoplanet ( Podolak 2003 , Movshovitz 
|et al.|[20T0} . We do not modify the gas opacities, but enforce a 
minimum for the grain opacities at the interface between the two 
regimes that corresponds to the gas minimum ensuring a smooth 
transition (see Ayliffe & Bate 2009a for more details). 



Table 1. Properties of the various models described in this paper, all of 
which were performed in a disc with a surface density of 750 g cnT 2 at r p . 
The opacity is given as a percentage of the interstellar grain opacity (IGO) 
assumed in the opacity tables we employ. Core radii are all based on the 
models of Seager et al. 2007 multiplied by factors of 3 and 10 as marked. 



2.5 Planetary Cores 

The planetary cores in these simulations are modelled by a grav- 
itational potential, and a surface potential that yields an opposing 
force upon gas within one core radius of the core's surface. The 
combination of the gravitational and surface forces takes the form 
of a modification to the usual gravitational force as 



r 2 



2R C -r\ 4 
R c 



(1) 



for r < 2 R c where r is the radius from the centre of the plane- 
tary core, R c is the radius of the core, and M c is the mass of the 
core (see Ayliffe & Bate 2009a for further details). This equation 
yields zero net force between a particle and the planetary core at 
the surface radius R c , whilst inside of the core's radius the force is 
outwards and increases rapidly with decreasing radius. Gas parti- 
cles therefore come to rest very close to the core radius, though the 
equilibrium position is slightly inward of this value due to the pres- 
sure exerted by the gas that accumulates on top of the inner most 
layer of particles. 

Seager et al. (2007) calculate core radii for solid exoplanets, 
amongst which are cases comprising of 75 per cent water, 22 per 
cent silicates, and 3 per cent iron. We use these models to deter- 
mine the realistic sizes of protoplanetary cores that correspond to 
the masses used in this paper. These core radii were employed for 
a number of calculations in Ayliffe & Bate (2009a), but these cases 
were not evolved for many orbits due to the short timesteps required 
deep within the planetary potential. To follow the evolution over 
longer periods it has been necessary to perform calculations using 
larger core radii. To this end we scale up these previously used core 
radii by factors of 3 or 10 to reduce the required computation time. 
We also performed new models using the realistic core radii, but 
these calculations were too slow to give useful results and thus are 
not reported here. The properties of our different models are given 
in Table [T] The rate at which the planetary core accretes gas may 
be affected by the form of the surface potential described by equa- 
tion[T| and this is explored briefly in Section |3~ij A smooth start to 
the calculations is provided by shrinking exponentially towards the 
desired R,. from an initial radius of 0.01 r p over the course of the 
first orbit. 
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Figure 1. Accretion histories for protoplanets with initial core masses of 
15 (A), 20 (D, E), and 33 ( J) M®. The models possess core radii 10 times 
their likely physical size. Model E is performed using a disc with 1 per cent 
interstellar grain opacity, whilst all the others use 0.1 per cent; model E 
is otherwise identical to D. The size and opacity dependencies are further 
illustrated in Fig. [2] 



2.6 Measuring the gas accretion rates onto the planetary 
cores 

We measure the gas accretion rates by calculating the rate at which 
gas passes into the self-consistently calculated Hill sphere of the 
protoplanet given by 



3M,' 



(2) 



where M p is the protoplanet mass which is the sum of the core mass 
CM C ) and the accreted mass (M acc ), where accreted mass comprises 
all the gas within R H . The gas mass is discretised amongst the SPH 
particles, allowing iteration through equation [2] until such time as 
the addition of a particle's mass to M acc no longer increases Ru 
sufficiently to encompass the next available particle. 

The net flux through the Hill radius corresponds to the growth 
rate of the envelope and/or circumplanetary disc, which are the only 
repositories for gas that fails to remerge from this region. Our use 
of the Hill radius to measure the mass growth is arbitrary, but rea- 
sonable, since any protoplanet must be smaller than the Hill radius. 
It is important to note that the Hill radius does not define the extent 
of the envelope, which tends to be smaller (e.g. ~ 0.25 R H Lissauer 
|et al.| |2009). However, the difference between the mass accretion 
rates (and total accreted masses) as measured at the Hill radius and 
at 0.3 Rh is small except at the start of the calculation (Section[3~2|l. 



3 RESULTS 

We have performed three-dimensional radiation hydrodynamics 
models of the accretion of gas onto planetary cores (or embryos) 
with a range of masses, over hundreds of orbits, in discs of varying 
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Figure 2. Accretion history for a 20 M ffi core with 0. 1 per cent and 1 per cent 
opacities (as marked), for two different core radii, 3x physical (dashed lines, 
models B & C), and lOx physical (solid lines, models D & E). Note that for 
the higher opacity case, which accretes more slowly, the models of different 
assumed core radii coevolve. For the lower opacity cases, which accrete 
large envelopes much more rapidly, the case with the smaller core radius 
accretes more slowly that the large core case during the initial ~ 50 orbits. 
However, thereafter, when their accretion rates are measured at equivalent 
masses, they exhibit very similar rates. 



opacity. This work extends upon models we performed in Ayliffe & 
|Bate| ( |20"09a) where we followed the accretion for a relatively short 
period. Moreover, by using discs with highly reduced opacities, the 
models presented here include the accretion of much more signif- 
icant envelope masses. In one case the accreted mass is sufficient 
to trigger a hydrodynamic collapse of the envelope (Section \3.2\ , 
followed by a return to steady gas accretion. 



3.1 Envelope accretion 

We performed calculations starting with various core masses rang- 
ing from 15 - 33 M ffi (Table[T}. The principal property that controls 
the accretion rate of a protoplanet of a given mass is the opacity of 
its envelope. A lower opacity enables more rapid radiative cooling, 
allowing the envelope to contract more quickly and so accrete gas 
at a faster rate ijH ubickyj et al. 2005 , Papaloizou & Nelso n|2005| 
Ayliffe & Bate 2009a I. The models presented in this work exploit 
this dependence to accelerate the growth process by adopting re- 
duced opacities. The impact of opacity on the growth rate of a pro- 
toplanet is demonstrated by comparing Models E and D in Figs.[T| 
and|2] Model E employs an opacity ten times larger than D, which 
results in a much slower rate of accretion for the former. 

In each case we introduced a bare core with no preexisting 
envelope. This results in a rapid initial gas accretion rate to form 
a quasi-static envelope, followed by a period of slowing accre- 
tion. This initial growth phase is not realistic for a protoplanet 
that forms in situ and concurrently accretes both solids and gas 
during the core formation phase (e.g. Poll ack et al.|1996| |Alibert| 
|et al.| [2005l. In models that account for both the core and enve- 
lope growth, the gas accretion rate tends to have an extended pe- 
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Figure 3. A comparison between results from the lower surface density, 
shorter evolution calculations presented in Ayliffe & Bate 1 2009a I and the 
longer, more rapidly accreting models discussed in this paper. The lower 
solid line is the accretion history over 160 orbits for a 33 M e (0.1 M^pi^) 
protoplanet in a disc with = 75 g cm" 2 , and interstellar grain opacity. 
Overlaid is a fit of the form M = aft + c (dashed line), where the fit is made 
over the range 10 - 50 orbits (shaded) and matches the measured growth 
over the entire evolutionary period. The upper line is the accretion history 
for a 33 M ffi protoplanet, in a disc of S p = 750 g cm" 2 , with 0.1 per cent 
IGO (model J), which by virtue of these conditions accretes much more 
rapidly. A fit to this curve is also shown. The departure of the evolution- 
ary model from the fitted curve is as discussed in Ayliffe & Bate 1 2009a I. 
The accretion rate initially declines as gas accretion builds a more optically 
thick envelope, and trapped heat prevents the envelope contracting, slowing 
its further accretion. However, when a sufficient mass is accreted, gravity 
comes to dominate the growth process, resulting in an accelerating accre- 
tion rate. 

riod during which it is almost constant. However, such an initially 
declining accretion rate was also seen in the semi-analytic models 
of Papaloizou & Nelson ( 2005 I, where the planetary accretion rates 
initial decreased due to the significant liberation of binding energy 
as mass fell into the planet's potential. This energy release heats the 
forming envelope, increasing the pressure, and reducing the rate 
at which mass can be added. Eventually, however, in all models, 
once a critical mass is reached, the energy release within the enve- 
lope is insufficient to further retard the accretion, which then accel- 
erates; this marks the transition from the thermally-dominated to 
gravitationally-dominated accretion regime. 

In |Ayliffe & Bate| < |2009a) , we studied the initially declin- 
ing accretion rate in one case which we followed for 160 orbits, 
but which did not reach the gravitationally-dominated regime. We 
found that the slowing envelope growth could be represented very 
well with an analytic fit of the form M = at + c where b = 0.40. 
In our new models, Model J is similar to the model from Ayliffe & 
|Bate|j2009a| l, but with a much lower opacity and a higher disc sur- 
face density. This produces higher accretion rates and, thus, allows 
us to reach the gravitationally dominated growth regime. We refit 
the old model using only orbits 10 - 50, and obtain essentially the 
same fit we obtained in Aylif fe & Bate| ( |2009a) using all 160 orbits; 



this is evident in Fig.|3]which illustrates that the fit is still extremely 
good at 160 orbits. Fitting Model J over the same 40 orbits, where 
an exponent b = 0.45 is found to be best, the trend of an initially 
reducing rate of accretion falters at the limit of the fitting region, 
around 50 orbits. At this point the accreted mass is a little more 
than half the core mass, and gravity begins to dominate the growth, 
leading to accelerating accretion. 

In all cases, such as the models shown in Fig.[JJ the transition 
from a slowing rate of accretion to an accelerating rate of accretion 
occurs once the envelope mass has reached half the core mass. This 
begins the process of runaway accretion, that continues until such 
time as the envelope potentially undergoes a dynamic collapse, as 
will be discussed in the following section. It is this period of accre- 
tion that leads to the most significant growth of the protoplanet's 
mass. 

As mentioned in Section [23| in order to make these three- 
dimensional models computationally viable, we were forced to 
adopt non-realistic core radii for the planetary surfaces. In Ayliffe 
& Bate (2009a I we compared the accretion rates achieved in mod- 
els using different core radii and found that for high opacities the 
accretion rates did not depend significantly on the core radius that 
was used, but for low opacities ( 1 per cent and 0. 1 per cent interstel- 
lar grain opacities) the accretion rates obtained with smaller cores 
were significantly lower than for equivalent models adopting larger 
cores. However, the earlier models were only evolved for 10 orbits, 
which is still during the initial phase of envelope creation when the 
accretion rate is decreasing. Fig.|2]illustrates the growth of 20 M e 
cores embedded in a disc of either 0. 1 or 1 per cent interstellar grain 
opacity (as marked). The solid lines present models with planetary 
cores 10 times the realistic core radii, and the dashed lines 3 times. 
Over the course of the initial settling period the different core radii 
in the 0.1 per cent opacity models lead to some divergence in the 
growth, as is evident. However, measuring the accretion rates of 
the two calculations at equivalent masses beyond the initial 50 or- 
bits, it is found that these rates deviate by no more than 20 per cent 
from one another. In the 1 per cent opacity case, the results ob- 
tained using the two different core radii are indistinguishable from 
each other. We believe that measuring the accretion rates at longer 
times, when the models are more established is responsible for the 
relatively small differences seen with different core radii. However 
in this case we are only varying the core radius by just over a factor 
of 3, whilst in Ayliffe & Bate (2009a I the comparison was made 
spanning more than a factor of 12. At the present time we do not 
have any models with which we can ascertain to what extent the 
accretion rate differences measured in this previous work were due 
to the large factor difference in the core radius, and what fraction 
came of the early times at which the measurements were made. 
As such, the enlarged radii of the cores that we have adopted here 
should be kept in mind. 

In the rest of this paper we focus on one particular case, Model 
J, that accretes a very significant envelope in a short period of 
time that undergoes a dynamic collapse. The other calculations dis- 
cussed in this section are ongoing, and will eventually allow us to 
explore the evolution of envelopes that are built up under less ex- 
treme conditions (i.e. with slower accretion rates). 

3.2 Hydrodynamical collapse - Model J 

Model J accretes the most significant mass of any of our models, as 
is to be expected given its favourable conditions. The 33 M e core is 
the most massive that we employ, and in this case is coupled with 
a large 10 times realistic core radius, and a 0.1 per cent interstellar 
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Figure 4. The accretion history of the gaseous envelope in model J, where 
the accreted mass is that contained within the self-consistently calculated 
Hill radius (solid line). Also shown is the mass evolution with time within 
radii of 0.3 (dashed), 0.1 (dot-dash), and 0.03 Rh (dots-dash). The envelope 
undergoes a hydrodynamical collapse after around 190 orbits when its mass 
is as 0.375 Mj u pj ler , with the resulting shock pushing material out of the Hill 
radius, whilst the remaining mass becomes more centrally condensed. This 
phase is shown in the inset panel, which illustrates the central condensation 
by the increasing mass within small radii as the overall mass falls. Follow- 
ing the collapse, accretion resumes, replacing the mass lost due to the shock 
propagation. 



grain opacity. The mass evolution of this protoplanet is shown in 
Fig. [4] As in all cases, the accretion rate initially drops away over 
the first ~ 50 orbits as the model settles, before beginning to in- 
crease. The acceleration of accretion becomes more marked at the 
crossover mass of a; 0.1 M Jupitcr , at around 100 orbits and continues 
during the runaway growth phase. A change occurs at around 190 
orbits, when the accreted mass has reached as 0.375 MjupiKr, At this 
point the mass of the envelope is such that it overcomes the pres- 
sure gradient which hitherto had supported it, and only by centrally 
condensing can a new hydrostatic equilibrium be established. To 
achieve this the envelope rapidly collapses, in the process trigger- 
ing a shock wave which propagates outwards from the very dense 
core that forms. This shock pushes material outwards, removing 
a small fraction of the material from within the Hill radius very 
rapidly at « 190.5 orbits. The drop in the mass within 1, 0.3, and 
0. 1 Rh can be seen clearly in the inset panel of Fig. [4] Most interest- 
ingly, the reduction of the mass within 0. 1 Rh indicates that there 
is a small loss of bound material as this radius falls well within the 
expected limit of a protoplanetary envelope. The magnitude of the 
mass loss from each volume is ~ 1 per cent of the mass therein. 
Following a brief increase (» 190.6 - 191 orbits), the mass within 
Rh continues to decline over the course of a further orbit before 
it begins to grow once again. However, at the same time the mass 
within 0.1 Rh increases rapidly, showing that the available material 
is quickly condensing around the core. Beyond ss 192 orbits the 
accretion rate measured at the Hill radius resumes its pre-collapse 
rate of 4 x 10~ 4 Mj up iter year . Moreover, this rate appears to be 



consistent at nearly all radii, indicating that gas is falling directly 
onto the new denser envelope, rather than becoming suspended in 
a more extended structure as was the case prior to collapse where 
the accretion rates were different at different radii. 

The process of envelope collapse is presented in Fig. [5] the 
left hand panels of which shows cross-sections in density in the 
Z-X plane through the centre of the planet, from a time preceding 
the collapse in the uppermost panel, and at various stages during the 
collapse in subsequent panels. In the first panel the density contours 
are near spherical at small radii, elliptical at 0.5 Rh, and pinched in 
towards the planet's poles at large radii. The second panel illustrates 
the strong shock propagating outwards from near the planet's core, 
which at this point has not altered the gas structure near the Hill 
radius (marked with a dotted line). However, it is possible to see 
that the collapse has been greater at the poles, deforming the previ- 
ously elliptical contours. The third panel demonstrates the contin- 
ued propagation of the shock, which is elongated along the vertical 
axis due to the lower density of material above and below the plane 
of the disc, which allows the shock to propagate more rapidly in this 
direction. Contours which had been pinched in at the planet's poles 
are forced outwards as the shock passes, but as can be seen, inside 
they are already resuming their pinched structure. The bottom panel 
illustrates the resettled state of the protoplanet and its surround- 
ings. The pinching at the poles, which in the top panel was found 
beyond 0.5 Rh, now extends down to ~ 0.1 Rh, and the surround- 
ing medium now forms a circumplanetary disc, having previously 
existed as an ellipsoidal envelope. The central density of the proto- 
planet, that is the gas density near the core, was »5x 10~ 4 g cm" 3 
in the top panel, and has increased to as 5.5 x 10~ 3 g cirT 3 in the 
structure shown in the bottom panel following the collapse. 

Returning to the formation of a circumplanetary disc, Fig. ^il- 
lustrates the midplane and vertical density structures about the pro- 
toplanet, emphasising the altered state of the envelope. The diver- 
gence from a spherically-symmetric density distribution, defined 
here as a difference of greater than 10 per cent between the mid- 
plane and vertical distributions, occurs at a radius 5 times smaller in 
the post collapse state than in the pre-collapse state. This illustrates 
that the envelope undergoes a more significant change in its struc- 
ture vertically, than it does in the plane of the disc, and is similar 
to the results seen in Ayliffe & Bate (2009a I. This was investigated 
further in Ayliff e & BateR 2009b) where it was found that circum- 
planetary discs formed around massive protoplanets, and that these 
discs tended to be thick, with dimensionless scale heights generally 
larger than 0.2. We measure the circumplanetary disc scale height 
by taking radial bins within the region with r < Rh/3 (measured 
from the protoplanet). In each radial bin, we take the SPH particle 
densities (from all azimuthal angles) and fit Gaussian profiles to the 
resulting vertical density distributions. A Levenberg-Marquardt al- 
gorithm is used to perform the fit, allowing the scale height to vary. 
This gives a measure of disc scale height versus radius, which in 
this case yields values for HJr that increase from 0.4 to 0.5 over 
the radial range of 0.05-0.33 Rh over which the disc extends. This 
radial extent of the circumplanetary disc can be seen in Fig. [7] in 
which the disc edge is taken to be at the location where the peak of 
the specific angular momentum occurs. Beyond this radius, the spe- 
cific angular momentum decreases with radius since, in the frame 
rotating with the protoplanet, the material in the circumstellar disc 
is counter-rotating relative to the gas captured by the protoplanet. 
Over the radial range 0.07-0.3 the specific angular momentum 
is on average 0.65 of the Keplerian values (marked by the dashed 
line). The displacement is due to the pressure support within the 
thick circumplanetary disc, and the degree of this displacement can 
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Figure 5. Cross-section plots illustrating the envelope collapse about a 33 M e (» 0.1Mj up i ler ) core. The Hill radius is marked as a dotted line. The accreted 
mass (gas within Ryj) is ss 0.375 Mj U pj tel prior to the collapse, reducing by x 1.3 per cent as the shock pushes material away from the core. The x and z axes 
are given in units of r p . Left panels: Cross sections in density illustrate the envelopes collapse, the shock propagation, and the formation of a circumplanetary 
disc. The central density increases by an order of magnitude from the top panel to the last. Right panels: Temperature cross-sections at equivalent times to the 
density panels. The peak temperature increases by 2100K to 6800K from the first to the last panel. The protoplanetary disc scaleheight is ss 0.05 r p , thus there 
is little material obstructing the vertical propagation of the shock front. As a result, the shock propagates more easily in the vertical direction than through the 
denser midplane, yielding the non-spherical propagation most clearly shown in the third righthand panel. 
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Figure 6. The density distribution along the x-axis (solid lines) and along 
the z-axis (dashed lines) through the protoplanet before (top panel) and af- 
ter (bottom panel) collapse. The plus symbols mark the radius and density 
at which the midplane and vertical density distributions differ by more than 
10 per cent. The vertical dotted line marks the protoplanet's Hill radius. The 
density distribution is modified significantly in both the midplane and ver- 
tically following collapse, but the point at which these distributions diverge 
moves inwards by a factor of 5 in radius. 



be used to approximately calculate the disc scale height as another 
check on the values measured above. Using the ratio of the specific 
angular momentum (J) to the Keplerian value (j k ) we obtain a scale 
height of 0.55. This is obtained using equation[3](see appendices B 
& C of Laibe et al. 2012), where we have used values of 1/2 and 
7/10 for the surface density and temperature exponents (p and q) 
for the circumplanetary disc, typical values from Ayliffe & Bate| 
(2009); the calculated scale height is not enormously sensitive to 
variations in these values within a reasonable range. 



H 



2(1 - jljk) 
p + q/2 + 3/2 



(3) 



This resulting scale height is somewhat larger than that which we 
directly measured, but has been calculated assuming a vertically 
isothermal disc and a lack of self-gravity, and is thus only approxi- 
mate. The specific angular momentum distribution further supports 
our assertion that a circumplanetary disc has formed as a result of 
the envelope collapse. 

The right hand panels of Fig. [5] shows the temperature struc- 
ture in a Z-X slice at equivalent times to those shown in the density 
panels to the left, allowing us to see the temperature evolution dur- 
ing the collapse and shock propagation. A hot front associated with 
the shock can be seen expanding away from the core with time, 
and the vertical elongation of the shock can be clearly seen in the 
third panel. Near the protoplanet's core, the peak temperature has 
increased by 2100K to more than 6800K in the post-collapse state 
shown in the final panel, and continues to increase for the remain- 
der of the calculation when accretion has recommenced. 

Fig.[8]shows the spherically averaged radial density and tem- 
perature profiles running outwards from the protoplanetary core at 
a number of points in the calculation. The change between 50 and 
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Figure 7. The specific angular momentum of the gas surrounding the proto- 
planet that comprises the circumplanetary disc. The vertical dot-dashed line 
marks Rh/3, the analytically expected edge of the disc, and this matches 
the measured turnover very well. The dashed line marks the Keplerian or- 
bital velocity based on the mass within the associated radius. Pressure sup- 
port within the disc means that a sub-Keplerian orbital velocity is expected, 
though with a similar gradient if the disc is rotating about the planet. This 
gradient matches reasonably between 0.07-0.3 Rh- 



150 orbits (the dashed lines) is relatively small, increasing as would 
be expected for a growing protoplanet. Arriving at the pre-collapse 
state (thin solid line) at around 190 orbits, at which point the accre- 
tion rate is at its maximum, the density and temperature maxima 
have increased more over 40 orbits than in the preceding 100 or- 
bits. Moreover, the temperature structure shows a marked change 
in its form. However the most significant changes occur during the 
collapse and 5 further orbits of evolution (thick solid line). The en- 
velope's collapse occurs very rapidly, and only stops when the new 
structure of the envelope is able to reestablish hydrostatic equilib- 
rium. For this to occur the density structure changes to deliver a 
much steeper gradient away from the planet's solid surface. This 
leads to a much steeper pressure gradient in this region, as shown 
in the bottom panel of Fig. [8] eventually satisfying the requirement 
VP = —pV<j>, where <p = GM(r)/r. As such, the envelope is able 
to resume steady accretion as the structure is able to bear the in- 
creasing weight; as mentioned previously, accretion resumes at the 
pre-collapse rate. 

Fig. [9] depicts the changing distribution of mass in the inner 
envelope from its pre-collapse state, to its post-collapse state. The 
total time span shown is 229 days, the time between pre and post 
collapse within the inner envelope, with these states marked using 
dashed lines. However, the most significant changes occur over less 
than 22 days at this scale, and this period is broken down in Fig. [5] 
into 4.3 day increments which are marked with solid lines. The Hill 
radius just prior to the collapse is equal to 0.054 r p , which taking 
|Lissauer et al.| p009 l's estimate of a ~ 0.25 Rh envelope radius, 
gives a size of 0.0135 r p . This radius is in reasonable agreement 
with the region over which the mass is significantly redistributed 
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Figure 8. Spherically averaged density (upper panel), temperature (middle 
panel), and pressure (lower panel) distributions about the protoplanet. The 
dashed lines are at times of 50, 100, and 150 orbits, their order ascending up 
the left hand axis. The solid lines show the pre-collapse state (thinner line), 
and the post-collapse state (thicker line). The pre-collapse state is equivalent 
to that in Fig. [9] whilst the post-collapse state is instead taken at the very end 
of the calculation in this case, when accretion has resumed its pre-collapse 
rate. 

during the collapse, which can be seen in Fig. [9] to be 0.01 r p (as 
0.2 R H ). 

Collapse of the protoplanetary envelope leads to a substantial 
increase in the temperature near the protoplanetary core, as dis- 
cussed above. A result of this temperature rise is an increase in the 
dissociation fraction of molecular hydrogen about the core, and an 
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Figure 9. Cumulative mass distribution calculated outwards from the core 
of Model J, where M acc = 0.375 Mister. The mass distribution is shown 
over times ranging from just before the collapse, to after the structure sta- 
bilises. The collapse proceeds very rapidly, as shown by the solid lines with 
cover a period of less than 22 days from first (lowest) to last (highest). 
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Figure 10. Fraction of atomic hydrogen versus radius within the inner re- 
gion of the protoplanetary envelope before it collapses (dashed line), and 
after it collapses (solid line). The higher temperatures that develop within 
the deep envelope when it collapses lead to a higher dissociation fraction of 
molecular hydrogen, whilst this process of dissociation will absorb energy, 
reducing the maximum temperature that is achieved. Prior to collapse the 
fraction of atomic hydrogen peaks at 0.33, whilst in the immediate after- 
math it is as high as 0.48. 
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enlargement of the region within which hydrogen is substantially 
dissociated; this can be seen in Fig.[l0] It is not the dissociation that 
triggers the collapse of the envelope, despite the process acting as 
an energy sink; the capacity of dissociation to absorb energy does 
lead to a lower final temperature within the envelope than might 
otherwise have been achieved. 



3.3 Accretion flow 

Whilst we find that a protoplanet envelope increases in mass 
throughout its evolution, excepting a brief period following a dy- 
namical collapse (seen at radii of 1, 0.3, and 0.1 Rh in Fig.|4j, it is 
not obvious that this accretion is a spherically symmetric process. 
Machi daet al.| ( [2008| > and |Tanigawa et aT7| ( |20 1 2[ > have performed 
three-dimensional calculations of protoplanet growth, and find that 
gas flows outwards along the midplane from a growing protoplanet, 
such that accreted material must be delivered vertically. In Ayliffe 
& Bate ( 2009b I we found that mass predominantly entered the Hill 
sphere along the midplane, but this analysis was simplistic in that 
it only considered inflow, failing to consider the possibility of out- 
flow, and we worked solely with the azimuthally integrated values. 
Here we make a more thorough assessment of the mass flow, and 
are able to look at this flow in a large extended envelope that has not 
undergone collapse, as well as in the dense envelope formed sub- 
sequent to such a collapse. It is the latter case which most closely 
resembles the principle model of Tanigawa et al. (2012) for a high 
mass protoplanet. 

Before and after the envelope collapse, we see material flow- 
ing in and out at the Hill radius in an alternating pattern, as can 
be seen in the first panels of Figs.[TT|&[T2] where negative values 
correspond to inflow, and the solid contour line denotes a flux of 
zero. This flow is easily explained as the passage of gas passing 
the planet and being deflected by the spiral shocks that form due to 
the gravitational perturbation provided by the protoplanet. Another 
contribution comes from gas following horseshoe orbits which also 
enter and leave the Hill sphere at broadly similar longitudes about 
the planet. Fig.[T3]ilfustrates the vector field that results in the mass 
flow observed at the Hill radius, which is marked in this figure with 
a dashed line. The data presented in Figs.^^&^Jwas constructed 
by considering the flow over a period of 4 orbits preceding and fol- 
lowing the collapse respectively. The stability of the gas flow about 
the protoplanet over these periods leads to the regular pattern that 
is seen at the Hill radius in both figures. In the pre-collapse case the 
pattern persists down to 0.3 as the second panel of Fig.pTjillus- 
trates. At smaller radii the flow is inwards across the spherical shell, 
with a larger flux at smaller radii as is expected for a consistent 
mass flowing across a decreased surface area. In the pre-collapse 
state the mass flow at small radii does not appear to possess any 
latitudinal dependence, rather it flows in almost spherically sym- 
metrically. Note that the figures are noisey at the poles as a result 
of the spherical polar grid used to calculate the flux, which leads to 
very small bins at high latitudes. 

In the post-collapse case, shown in Fig. [T2] the mass flux at 
the Hill sphere is little changed from the pre collapse case, and in 
both cases it is reminiscent of the structures seen in Fig. 5 of |Tani-| 
|gawa et aL 1 2012 1; note that we are plotting a time average flux, 



whilst 



Tanigawa et al. plot an instantaneous flux (pv r ). Unlike [Tani-| 



tude, with more significant outflow at longitudes of ~ 20 and 200 
degrees. This outflow originates at a radius of « 0.17 R H , which 
is around half the radius of the circumplanetary disc, inside of this 
radius the flow is inwards. Meanwhile at 0.3 ^h, mass is flowing 
inwards at higher latitudes. The relative fluxes are such that the net 
flux at each radius is negative, enabling the protoplanet to continue 
to grow. This growth is corroborated by Fig. [4] which shows the 
mass evolution of Model J within the 4 radii considered here, and 
indicates that the mass consistently increases within 0.3 Rh (dashed 
line) once the post-collapse state is achieved. Following the enve- 
lope's collapse, there is more discernible structure to the mass flow 
at the smaller radii. At 0.1 the inflow along the circumplanetary 
disc midplane is not as significant as the flow at high latitudes, as 
is shown in the third panel of Fig. [12] This panel marks a clear di- 
vide in the motion of gas at and around a radius of 0.1 ^h, where 
the flow is only found to be inwards, signalling that material within 
this radius is truly bound to the protoplanet. This is of particular 
interest because of the gas flow found at 0.03 Rh, and shown in 
the final panel of Fig. [T2] At this small radius there are signifi- 
cant flows of material, pushing outwards at intermediate latitudes 
(~ 45-70 degrees), and pouring inwards again along the midplane. 
This is indicative of significant circulation of the bound material 
below 0.1 Rh, and will be discussed further in Section|4~2"l 



4 DISCUSSION 



4.1 Hydrodynamic collapse 



|gawa et al.H2012| , who see these alternating structures persisting 
down to very small radii, our model has already lost any sign of 
the in-out flow pattern at a radius of 0.3 ^h- Instead, at 0.3 we 
find an outflow along the circumplanetary disc midplane at every 
longitude, though this outflow shows a slightly alternating magni- 



From the earliest suggestion of Perri & Cameron ( 1974) it has been 
thought that a giant planet might form through the hydrodynamic 
collapse of a gaseous envelope onto a solid core which caused it 
to assemble. This was followed by numerous models that effec- 
tively sought for hydrostatic solutions to various combinations of 
properties to establish when such a collapse might occur (jMizuno 
|et al.|1978[|Mizuno|1980"l|Sasaki|1989l ). The first models that at- 
tempted to model giant planet growth from the initial core forma- 
tion, through to the envelope growth were performed by |Boden"^] 
|heimer & Po llack ( 1986). These models revealed that a protoplan- 
etary envelope would gradually contract as the planet grew, leading 
to a quasi-static contraction beyond previously calculated values 
for the critical mass, as long as the envelope did not effectively de- 
tach from the protoplanetary disc (that is there was a sufficiently 
rapid supply of material from the latter to the former). Under these 
conditions, their was no evidence to suggest that the hydrostatic 
balance should reach some limit beyond which a collapse was in- 
evitable, and later semi-analytic models originating from these ear- 
lier works, suc h as|Pollack et aL| (1996 1, Hubickyj et al. (2005](, and 
|Lissauer et al.H2009| , suggest no need for a dynamic collapse. Our 
Model J follows a pattern of stable growth for the vast majority of 
its history, though not evidently undergoing any significant contrac- 
tion, and resumes this pattern of growth subsequent to its envelope 
collapse. 

The models presented in this article are performed using a 
three-dimensional hydrodynamics code that include self-gravity, 
and radiative transfer, but which omits the core formation phase, 
and the deposition of energy due to planetesimal accretion that are 
included in the semi-analytic works discussed. However, at the time 
of interest around the envelope collapse, the energy release is ut- 
terly dominated by the contraction of the gaseous envelope, such 
that solids accretion energy may be regarded as negligible. The 
metallicity of the envelope might be significantly modified by the 
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Figure 11. Mass flux through shells of various radii (marked in panels) surrounding the protoplanet core in Model J prior to its envelope collapse. At one Hill 
radius, the flow takes on a form that is similar to that seen in Fig. 5 of Tanigawa et al. 1 20121, and reflects the combined effect of horseshoe orbits and bent 
flow lines passing the protoplanet. Material both enters and leaves the sphere predominantly at the midplane where densities are highest; the marked contour 
line denotes a flux of zero such that the area within the contour marks the outflow. For ease of comparison with Tanigawa et al. 1 2012 1 the solar and anti-solar 
points are at longitudes of and 1 80 degrees respectively, and inflowing material is shown by a negative flux. At the smallest radii the flow is largely inwards 
across the shell, with an average flux of -1 X 10~ 2 code units (-5.5 X 10~ 5 g cm~ 2 s~'). 
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Figure 12. Identical to Fig.^^except illustrating the gas flow after Model J undergoes its envelope collapse. The flow at one Hill radius is very similar to the 
pre-collapse case, however at radii between 0.1 - 0.3 Ru the outflow is concentrated along the midplane, whilst inflow occurs at higher latitudes, reflecting the 
vertical accretion seen by Machida et al. 1 2008 1. At the smallest radii, there is evidence of a gas turn over, where material is flowing in along the midplane and 
at the poles, and out in a range of moderate to high latitudes (see Section[4!2). 
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Figure 13. Velocity vectors in the plane of the disc, illustrating the gas flow 
about the Hill radius (marked with a dashed line) that leads to the alternating 
pattern of in and out flow seen in the first panels of the mass flux plots 
shown in Figs. |ll| & |12| This vector field is plotted for a time preceding the 
envelope collapse, but a very similar field exists after the shock associated 
with the collapse has passed out of the region. 



ablation and evaporation of grains that have been accreted over the 
planets history, and we make no attempt to account for this. The 
opacity in Model J is reduced by a fixed factor of 10 3 , at the lowest 
end of the range suggested by Movsho vitz et al.| ( |2010] l, who found 
such opacities due to grain settling and coagulation in regions of 
the envelope. 

It is difficult to disentangle the causes and effects of the very 
rapid collapse we find, for example the surge in temperature leads 
to a higher fraction of dissociated hydrogen. However, as stated 
in Section [3~2| it does not appear that the dissociation of molec- 
ular hydrogen acts to trigger the collapse, as the fraction remains 
steady in the preceding period. There is also no evidence of the 
Kappa-mechanism acting within the protoplanetary atmosphere in 
our model, as was found by Wuchterl ( 1991 1 to cause a dynamic 
collapse. Wuchterl also found that this collapse led to a signifi- 
cant ejection of material from the protoplanet, whilst we find only a 
small drop of ss 1.3 per cent in the mass within the Hill radius, and 
a rapid increase at and below 0. 1 as the protoplanet structure 
shifts to its new state. 

It appears that our Model J protoplanet reaches the hydrostatic 
limit for its formative structure, and that the internal pressure gra- 
dient can no longer accommodate the addition of mass by a small 
adjustment. This is illustrated in Fig. [14] which shows the ratio 
of the pressure force to the gravitational force against radius at a 
number of stages of the envelope collapse. From an initial state, 
in which the protoplanetary atmosphere is in hydrostatic equilib- 
rium out to a radius of 0.0067 r p (0.12 Ru), the atmosphere rapidly 
begins to restructure, pushing the hydrostatic region down to a ra- 
dius of ss 0.0013 r p (0.025 Ru). This initial stage leaves the form 
of the graph otherwise relatively unchanged, but as the central con- 
centration of mass continues, a shock begins to form as the pres- 
sure near the core surges. The maximum pressure within 0.001 r p 
has increased by an order of magnitude between the first panel and 
the third, leading to a somewhat steeper gradient over this region. 
However, at this point in time the gradient between 0.001-0.002 r p 
steepens much more rapidly, forming a shock front, and this front 
marks the new radial limit of hydrostatic equilibrium as can be seen 
in the third panel. The subsequent panel shows the shock prop- 
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Figure 14. Each panel shows the ratio of the pressure force to the gravita- 
tional force in the envelope and beyond, with the first and last panels corre- 
sponding in time to the pre and post-collapse cumulative mass distributions 
shown in Fig. [9] A ratio of one indicates that the material is in hydrostatic 
equilibrium, which before the envelope collapse applies to a region out to 
0.0067 r p (0.12 Ru, top panel, marked with vertical dotted line), but which 
post-collapse reaches out to just 0.0009 r p (0.016 Ru). By the final panel 
the shock wave has cleared the inner 0.01 r p , leaving the environment in- 
ternal to this in a new hydrostatic equilibrium, whilst its mass continues to 
increase, as shown in Fig. [4] 
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agating outwards, whilst the hydrostatic core shrinks a little. By 
the final panel the inner envelope has resettled and the structure 
has stabilised, whilst the shock's propagation continues outwards, 
eventually moving beyond the limits of the modelled region. 

It is possible that this envelope collapse only occurs due to the 
high accretion rates achieved due to our selected disc conditions. 
The low opacity assumed promotes very rapid planet growth, and 
it may be this rapidity that prevents the envelope from adjusting 
its structure more gradually to accommodate the increasing mass. 
As such, it may be that the hydrodynamical collapse of a plan- 
etary atmosphere can only occur if that planet if accreting very 
rapidly. Further models will be required to determine whether or 
not this is the case. We note however that the accretion rate mea- 
sured in Model J, both just before and after the collapse, is still 
not as rapid as would be found using a locally-isothermal equa- 
tion of state, despite the large reduction in opacity. Lissau eFet al.| 
(2009 ) present accretion rates in their fig. 3, where these rates were 
obtained by |D'Angelo et al.| {2003b[ l using three-dimensional hy- 
drodynamical models with a locally-isothermal equation of state. 
Applying our disc conditions to their results yields an accretion 
rate of 8 x 10~ 4 M Jupiter year -1 , which is twice the rate measured in 
our radiation hydrodynamics models prior to collapse. 

At this juncture we note that the results given in |Lissauer| 
|et al.| ( |2009| > show a viscosity dependence, where higher viscosi- 
ties lead to more rapid accretion. In our calculations viscosity is 
not a constant, but is proportional to the spatial resolution of the 
SPH method. Thus, the viscosity is lower in regions of higher den- 
sity, and these differences mean the above comparison is only ap- 
proximate. Further, we note that in the absence of a protoplanet, 
the unperturbed circumstellar disc in our models has a viscosity of 
a x 4 x 10~ 3 , consistent with the fixed viscosity global models of 
|Bate et aI7| < |2003^ that are used to inject gas at the boundaries of 
the disc section. It is these boundaries that determine the rate at 
which gas is supplied to the disc section in these local models. As 
such, once the disc is perturbed, the spatially varying viscosity of 
the SPH calculations leads to to an inconsistency with the global 
models, and so the boundaries. A further caveat arising from the 
boundary implementation is that the injected material comprises 
gas on both circulating and librating orbits (Lubow et al. 1999J, or- 
bits that are modified as the gas passes through the modelled disc 
section. However, these modifications are lost when the gas leaves 
the section, and new gas is injected without these modifications, 
leading to a further inconsistency. 

4.2 Atmospheric turn over 

As briefly mentioned in Section|33]in reference to the third panel 
of Fig. ^] there appears to be significant motion of the bound gas 
besides rotation in the plane of the disc. An apparent rolling mo- 
tion is indicated by the flow through the surface at 0.03 Ru, the 
gas unable to escape, but flowing out and in within the bound at- 
mosphere. These flows appear to be convection cells in the deepest 
parts of the planetary envelope where the medium is extremely op- 
tically thick, and thus unable to cool readily by radiation. Fig.|15| 
shows a region out to 0.07 R H (3.7 x 10~ 3 r p ), at which radius the 
temperature is 675K, increasing to a maximum of 7400K near the 
core, demonstrating a significant thermal gradient against radius. 
Bodenheimer & Pollack ( 1986) found strata of convection within 
protoplanetary envelopes in their one-dimensional models, with a 
dense inner convection region that contained some 95 per cent of 
the envelope mass, and depending on the age and core size, mul- 
tiple higher level convection zones were also found to form. How- 
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Figure 15. A slice through Model J after the envelope has collapsed, show- 
ing the mean velocity vector field over the course of 4 orbits which corre- 
spond to those used to produce the mass flux plots in Fig. [T3] The outer 
circle corresponds to the spherical surface of Fig.rT^jat 0.03 Ru, whilst the 
inner filled circle illustrates the radius of the planetary core. This slice in 
the z-r cy i plane rotates about the core at the same average angular velocity 
as the gas, which exhibits solid-body rotation within 0.04 Ru, so that the 
vector field is seen in the rotating frame of the gas. The bold sections of the 
outer circle cover the latitudes of 45 - 70 degrees at which out flow was 
seen in the final panel of Fig. |12| 



ever, it is likely that the convection in fully 3-D models that include 
rotation, is very different. 

Fig.^Jillustrates the vector field in a slice through the proto- 
planet, averaged over 4 orbits that correspond to those used to pro- 
duce the mass flux rendering of Fig.[T2] Within ~ 0.04 R H the proto- 
planet rotates as a solid body, with a rate of 4.65 x 10~ 7 rad s -1 (giv- 
ing a day equivalent to as 160 Earth days). The plane to which the 
velocity vectors have been mapped was rotated at this rate, such that 
the inner region of the atmosphere remained consistently aligned 
with the plane over the time of averaging. There are four distinct 
regions of turn over in which gas flows outwards at mid-latitudes 
before falling inwards again along the disc midplane. These looping 
flows are subsonic, with mean velocities of less than 3 m s -1 and a 
maximum of less than 8 ms" 1 in a region where the sound speed is 
greater than 6 km s -1 . A time average to a fixed plane, or a spatial 
average rotating about the z-axis at a single moment in time both 
reveal a similar structure, suggesting that the general form of these 
structures is persistent. The bold sections of the 0.03 Ru circle are 
marked between 45 - 70 degrees, the latitudes where outflow was 
seen in the fourth panel of Fig.[T2] As might be expected, the veloc- 
ity vectors across these bold segments indicate outflow, illustrating 
that these looping flows are responsible for the features in the ren- 
dered plot. The 0.03 Ru region marked by the outer circle contains 
« 97 per cent of the total protoplanetary mass, that is the mass mea- 
sured out to the Hill radius; prior to collapse the same region had 
contained just under 40 per cent of the mass within the Hill radius. 
In this dense environment, possessing a significant radial tempera- 
ture gradient, the 4 distinct cells revealed in the velocity field might 
well be indications of convection. 
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Figure 16. Longitudinally (azimuthally) integrated mass flux for Model J at the same four radii considered in Figs. |ll| & |12| The solid lines represent the 
pre-collapse case (Fig. |l 1) and the dashed lines the post-collapse case (Fig. 1 12) . In the former, it is notable that the only instance in which an outward flux 
can be seen is for the largest radii, where this is easily understood by consideration of the large scale flows (spiral shocks, horseshoe regions, and vertical 
accretion) surrounding the planet, and persists in a similar form post-collapse. Post-collapse there is more structure visible in the flow at small radii. At 0.3 Rn 
there is a midplane outflow that corresponds to that visible the corresponding panel of Fig. [12] whilst at 0.1 Rn all material is flowing inwards, though at the 
midplane the flux is small relative to other latitudes. Finally, at 0.03 i?n there is evidence of significant turnover in the bound envelope post-collapse, possibly 
indicative of convection. 



4.3 Gas accretion 

Tanig awa et al.| (2012 ) have recently examined gas flow onto and 
within circumplanetary discs. They found, in agreement with the 
work of Machida et al. (2008), that gas flowed onto a circumplan- 
etary disc predominantly in the vertical direction. Moreover, their 
works suggests that material is flowing out along the midplane of 
a circumplanetary disc. Once the envelope collapses in Model J of 
this work, and a circumplanetary disc forms, we also find that ma- 
terial is flowing outwards along the midplane beyond a radius of 
0.17 Ru, which is at around half the outer radius of the circum- 
planetary disc jQuillen & Trilling|| 19981 |Ayliffe & Bate]|2009b| 



Martin & Lubow|201 1|>. Fig. 16 which is equivalent to Fig. 6 in 



Tanigawa et al.| q2012 ), illustrates the longitudinally (equivalently, 



azimuthally) integrated mass flux for the pre (solid line) and post 
(dashed line) collapse flows shown in Figs. [TT]& [T2] Before the 
envelope collapses, at 0.3 Rn the outflow seen in two places at the 
midplane is counterbalanced by the associated midplane inflows 
and the inflow at higher latitudes, yielding a net inflow, as shown by 
the solid line in the second panel of Fig. [16] However, post collapse 
the consistent midplane outflow seen for this radius at all longitudes 
results in a peak of outflowing material in a region ±30 degrees lat- 
itude. 



About the planet, at one Hill radius the flow is dominated by 
the streams of material passing in and out of the region due to the 
form of their circumstellar orbits (see Fig. |13) . As such, it is un- 
surprising that the fluxes we obtain, and those of Tanigawa et al. 
( |2012fr scaled to our model, are similar at this radius; note we make 
the following comparisons using our post-collapse case which bet- 
ter resembles Tanigawa et al. 's models. At 1 R a they obtain a peak 
outflow along the midplane of 1.8 x 10~ 7 gcirT 2 s , whilst we 
obtain a value of 2.5 x 10~ 7 g cirT 2 s . The form is also similar, 
with wings of inflow at higher latitudes of similar magnitude; at 
the highest latitudes, our normalisation by area leads to a non-zero 
inflow. However, at smaller radii the mass fluxes we find are con- 
siderably larger than at the Hill radius, which differs substantially 
from |Tanigawa et aT] s results, where the peak flux at all radii dif- 
fer by less than an a factor of 4. At 0.3 Rn the peak outward flux 
has grown to 2.2 x 10 -6 gem -2 s -1 , and at 0.1 Rn the outflow has 
ceased, but the inflow flux at high latitudes has increased by an- 
other factor of around 5. The final panel of Fig.[l6]reveals the mass 
flow resulting from the formation of convection cells in the deep 
atmosphere post collapse. 
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5 SUMMARY 

We have performed three-dimensional self-gravitating radiation hy- 
drodynamical models of planet growth that, subsequent to an ex- 
tended period of growth, result in a dynamical collapse. A series 
of models were performed using different core masses, core radii, 
and opacities, that extend the range of accreted mass achieved in 
Ayliffe & Bate (2009a). We present the first results from a three- 
dimensional hydrodynamical model of planet growth by core accre- 
tion that has been found to produce a hydrodynamic collapse. The 
result of this collapse is a very centrally-condensed protoplanet, 
surrounded by a circumplanetary disc, that continues to accrete. 
The inner reaches of the envelope have undergone significant dis- 
sociation of molecular hydrogen, and appear to possess convection- 
like cells of gas turn over, whilst the inflow of new gas occurs near 
vertically at high latitudes. 

The circumplanetary disc, with radius » and dimension- 
less scaleheight of 0.4 - 0.5, exhibits a reversal in the direction of 
mass flow along the midplane at around 50 per cent of its radius; 
that is to say there is inflow only within the inner 0.17 R H . The 
degree of central condensation in the post-collapse state leads the 
model to show good agreement with previous calculations consid- 
ering mass flow that have presumed a pre-existing high mass core 
of relatively small size (Tanigawa et al. 2012 ). Conversely, before 
the collapse, the extended protoplanetary envelope exhibits a more 
spherically symmetric inflow of material. 

To achieve rapid growth in these models we have adopted very 
favourable conditions, particularly a low opacity of just 0. 1 per cent 
of the interstellar grain opacity. It may be that the hydrodynamic 
collapse found in this work is a result of the very rapid growth of 
the envelope, promoted by these disc conditions, though this can- 
not be said definitively without performing further models in less 
favourable discs. 
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